Stochastic Global Optimization =============================== Why stochastic global optimization ----------------------------------- DSGE and VAR likelihoods are routinely multimodal, with shallow ridges, flat plateaus, and parameter-region dependence on regime probabilities. Local optimizers (``fmincon``, ``csminwel``, ``newrat``) converge to the basin they happen to start in. Stochastic global optimizers trade per- iteration cost for coverage: they sample the bounded box ``[lb, ub]`` with a population of candidate solutions, exchange information across candidates, and contract toward the best region only as the iteration budget is consumed. They are the right tool when: - the model has more than a handful of parameters; - the user does not have a strong prior on a starting region; - a posterior maximization is intended to seed an MCMC sampler that must explore beyond a single mode. RISE ships a unified family of nine global optimizers (plus the legacy ``@bee`` and the deterministic ``mads`` baseline), all sharing the same calling convention and outer-loop machinery. The unified shell ------------------ Every algorithm in the ``+globalopt`` package is a thin entry point that wires an algorithm-specific *step function* into ``globalopt.generic_loop``. The shell handles: - initial-population construction; - iteration counting and budget checks (``MaxIter``, ``MaxFunEvals``, ``MaxTime``); - progress display (``Display = 'iter' | 'on' | 'off'``); - manual stopping via a ``ManualStoppingFile.txt`` sentinel in the current directory (see ``utils.optim.manual_stopping``); - optional finite-difference Hessian at the optimum (returned as the 4th output). This keeps every algorithm's implementation small: the step function implements only the algorithm-specific update; everything else is shared. Common calling convention ------------------------- All nine algorithms expose the same signature, both via the package and via the ``rise_`` shortcut:: [xfinal, ffinal, exitflag] = rise_lshade(fh, x0, lb, ub) [xfinal, ffinal, exitflag] = rise_lshade(fh, x0, lb, ub, options) [xfinal, ffinal, exitflag] = rise_lshade(fh, x0, lb, ub, options, hyper) [xfinal, ffinal, exitflag, H] = rise_lshade(...) Inputs: - ``fh`` -- function handle, ``fval = fh(x)``. - ``x0`` -- initial guess column vector. Used by some algorithms to seed the initial population (e.g. CMA-ES's initial mean); ignored by others. - ``lb``, ``ub`` -- column vectors with lower and upper bounds. - ``options`` (optional struct) -- the MATLAB-style fields below. - ``hyper`` (optional struct) -- algorithm-specific hyperparameters (see each algorithm's entry). Outputs: - ``xfinal`` -- best parameter vector found. - ``ffinal`` -- best objective value. - ``exitflag`` -- 1 on success. - ``H`` (optional 4th output) -- finite-difference Hessian at ``xfinal`` (computed lazily on request). Common options (every algorithm) ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ================ =========================== ==================== field meaning default ================ =========================== ==================== ``MaxIter`` outer-loop iteration cap ``200`` ``MaxFunEvals`` total function-eval cap ``inf`` ``MaxTime`` wall-clock cap (seconds) ``inf`` ``Display`` ``'iter'`` / ``'on'`` / ``'off'`` ``'off'`` ================ =========================== ==================== The algorithm catalogue ----------------------- The nine algorithms group naturally into four families. Each entry gives the philosophy, the seminal reference, and the algorithm-specific ``hyper`` fields. Differential Evolution family ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ These three algorithms share the DE skeleton -- a population of candidates, mutation by weighted vector differences, binomial crossover, greedy selection -- and add successive refinements. **DE** (``rise_de``, ``+globalopt.de``) The original Storn-Price algorithm (1997). Strategy selectable via ``hyper.strategy``: ``'rand/1'``, ``'best/1'``, ``'current-to-best/1'``, ``'rand/2'``, ``'best/2'``. Hyperparameters: ``NP`` (population, default ``20 + 5*dim``), ``F`` (mutation factor, ``0.8``), ``CR`` (crossover probability, ``0.9``), ``strategy`` (``'rand/1'``). **L-SHADE** (``rise_lshade``, ``+globalopt.lshade``) Tanabe & Fukunaga (2014), CEC 2014 winner. Three additions: JADE-style ``current-to-pbest/1/bin`` mutation with an external archive; per-individual ``F_i`` (Cauchy) and ``CR_i`` (Normal) sampled from a success-history memory (SHADE); linear population-size reduction from ``NP_init`` to ``NP_min`` as the budget is consumed. Hyperparameters: ``NP_init`` (``min(18*dim, 50)``), ``NP_min`` (``4``), ``H`` (history memory size, ``6``), ``p_best`` (greediness, ``0.11``), ``arc_rate`` (archive size, ``2.6``). **jSO** (``rise_jso``, ``+globalopt.jso``) Brest, Maucec & Boskovic (2017), CEC 2017 winner. Refinement of L-SHADE: weighted mutation ``F_w * (pbest - x) + F * (r1 - r2)`` with ``F_w`` depending on the fraction of the budget consumed (suppresses overly large early steps); ``F`` capped at 0.7 before 60% of the budget; generation-dependent greediness; smaller archive (``arc_rate = 1.0``). Hyperparameters: ``NP_init`` (``min(round(25*log(dim)*sqrt(dim)), 50)``), ``NP_min`` (``4``), ``H`` (``5``), ``p_max`` (``0.25``). Population / swarm metaheuristics ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ **ABC** (``rise_abc``, ``+globalopt.abc``) Artificial Bee Colony (Karaboga 2005, with the Akay-Karaboga 2012 per-dimension modification rate). Three roles per cycle: *employed* bees do one mutation per food source; *onlookers* run ``NP`` roulette-wheel mutations weighted by fitness; *scouts* uniformly replace food sources whose trial counter exceeds ``limit``. Hyperparameters: ``NP`` (food sources, ``max(20, 5*dim)``), ``limit`` (abandonment threshold, ``NP*dim``), ``MR`` (per-dimension modification rate, ``0.5``), ``phi_low``, ``phi_high`` (perturbation bounds, ``[-1, 1]``). **ACO** (``rise_aco``, ``+globalopt.aco``) Continuous Ant Colony Optimization, ACO_R (Socha & Dorigo 2008). Keeps a ranked archive of ``k`` solutions with Gaussian-kernel weights. Each iteration samples ``m`` new ants by picking a guide from the archive (roulette), sampling each coordinate from ``N(guide_d, sigma_d)`` with ``sigma_d`` set by the mean absolute distance in dimension d to the other archived solutions, then merging into the archive. Hyperparameters: ``k`` (archive, ``max(20, 5*dim)``), ``m`` (ants per iteration, ``k``), ``q`` (elitism, ``0.1``), ``xi`` (convergence, ``0.85``). **BBO** (``rise_bbo``, ``+globalopt.bbo``) Biogeography-Based Optimization (Simon 2008). Habitats are sorted best-first; the top ``n_elite`` are preserved; each remaining habitat receives features (decision-variable values) by migration from randomly chosen "donor" habitats with rates set by ``I`` (immigration) and ``E`` (emigration). Features are then perturbed by Gaussian creep with std ``mut_sigma`` (as a fraction of ``ub - lb``). Hyperparameters: ``NP`` (``max(20, 5*dim)``), ``I`` (``1``), ``E`` (``1``), ``mut_prob`` (``0.05``), ``mut_sigma`` (``0.02``), ``n_elite`` (``2``). Knowledge-based ~~~~~~~~~~~~~~~ **AGSK** (``rise_agsk``, ``+globalopt.agsk``) Adaptive Gaining-Sharing Knowledge (Mohamed, Hadi & Mohamed 2020, CEC 2020 winner). At every generation the population is sorted by fitness; for each individual a per-dimension partition splits the n coordinates into a *junior* subset (local exploitation using rank-neighbours) and a *senior* subset (global exploration using top/middle/bottom slices). The junior fraction shrinks as the budget is consumed. Update factors ``K_F``, ``K_R`` are sampled per individual from SHADE-style success-history memories. Hyperparameters are mostly self-tuning; the implementation here is the single-memory variant. Evolution strategy ~~~~~~~~~~~~~~~~~~ **CMA-ES** (``rise_cma_es``, ``+globalopt.cma_es``) Covariance Matrix Adaptation Evolution Strategy (Hansen 2016). The distinct one in the family: model-based rather than mutation-based. Maintains a multivariate Gaussian sampling distribution; updates its mean by weighted recombination of the best offspring, and its covariance via rank-1 and rank-mu updates that adapt to the local geometry. The default adaptation constants are derived from ``lambda`` and ``dim``. Hyperparameters: ``sigma0`` (initial step size, ``0.3``), ``lambda`` (population, ``4 + floor(3*log(dim))``), ``mu`` (parents, ``floor(lambda/2)``). Direct-search baseline ~~~~~~~~~~~~~~~~~~~~~~ **MADS** (``rise_mads``, ``+globalopt.mads``) Mesh Adaptive Direct Search (Audet & Dennis 2006), in the coordinate-aligned variant closer to Generalized Pattern Search (GPS) than to OrthoMADS. **Deterministic** given ``x0`` and the bounds: no internal randomness. Coordinate-aligned poll directions (``+/- e_i``), opportunistic polling (accept the first improving direction), per-dimension step sizes that contract on failure (``contract = 0.5``) and expand on success (``expand = 2.0``). Useful as a non-stochastic reference when comparing the population-based optimizers. Hyperparameters: ``init_step`` (``0.1*(ub - lb)``), ``tol_step`` (``1e-8``), ``expand`` (``2.0``), ``contract`` (``0.5``), ``opportunistic`` (``true``). Calling from ``estimate`` ------------------------- The ``'optimizer'`` option of ``estimate(...)`` accepts a function handle, a string, or a cell ``{name, options...}``. All ``rise_`` shortcuts are valid drop-ins:: m = estimate(m, 'data', mydata, 'optimizer', @rise_lshade); m = estimate(m, 'data', mydata, ... 'optimizer', {'rise_lshade', 'MaxFunEvals', 5e4}); To override the algorithm hyperparameters at estimation time, write a small wrapper that closes over the desired ``hyper``:: my_lshade = @(fh, x0, lb, ub, opts, varargin) ... rise_lshade(fh, x0, lb, ub, opts, struct('NP_init', 100)); m = estimate(m, 'optimizer', my_lshade); Calling directly ---------------- Outside ``estimate``, the algorithms are standalone optimizers and can solve any objective:: rosenbrock = @(x) sum(100*(x(2:end) - x(1:end-1).^2).^2 ... + (1 - x(1:end-1)).^2); n = 10; lb = -5*ones(n,1); ub = 5*ones(n,1); x0 = lb + (ub - lb).*rand(n,1); options = struct('MaxFunEvals', 1e4, 'Display', 'iter'); [xstar, fstar] = rise_lshade(rosenbrock, x0, lb, ub, options); Manual stopping --------------- Long runs can be stopped cleanly without losing the current best: create a file named ``ManualStoppingFile.txt`` (any content) in the working directory. ``utils.optim.manual_stopping`` is polled inside ``generic_loop``; the next iteration boundary terminates the run and returns the current best. The legacy ``@bee`` class ------------------------- RISE also retains the older ``@bee`` class (``m/optimizers/@bee/bee.m`` with ``m/optimizers/bee_gate.m`` as the function-style entry point). It exposes a different option set and remained the default for RISE's MCMC posterior-mode search until the ``+globalopt`` family was introduced. New work should prefer one of the nine ``rise_`` shortcuts above; ``@bee`` is kept for backward compatibility with existing scripts. Algorithms not currently included --------------------------------- Two algorithms sometimes asked about -- *firefly* (Yang 2008) and *glowworm swarm optimization* (Krishnanand & Ghose 2009) -- are not in the ``+globalopt`` family. The family is opinionated toward algorithms that have been independently winners or runners-up of the CEC real-parameter optimization competition (DE, L-SHADE, jSO, AGSK, CMA-ES), plus the most-cited swarm metaheuristics (ABC, ACO, BBO), plus a deterministic direct-search baseline (MADS). Firefly and glowworm slot in cleanly if needed: drop a ``+globalopt/+firefly/`` package with ``firefly.m`` (thin entry point wiring ``hyper.algorithm_step``) and ``firefly_step.m`` (init / iterate modes), matching the existing 9 algorithms' shape.